### CODE FOR:
###
### Formative Personal Experiences:
### How Benefit Recipiency and Income Changes Shape Perceptions of System Abuse
###
### VERSION: 2024-FEB-19


rm(list=ls())
gc()



# WORKING DIRECTORY -------------------------------------------------------
setwd("C:/Users/miros/Desktop/Research/Formative personal experiences")



# LIBRARIES ---------------------------------------------------------------
library(haven)
library(fixest)
library(tidyverse)
library(tableHTML)
library(RColorBrewer)
library(ggtext)
library(summarytools)
library(ggpubr)
library(cowplot)




# DATA --------------------------------------------------------------------

# __ Loading --------------------------------------------------------------
# Note: Data are available for research purposes, including replication, 
#       on the webpage of Norwegian Centre for Research Data (SIKT) at 
#       https://doi.org/10.18712/nsd-nsd2458-v6

SUPPA <- subset(read_dta("02-data/SuppA_1-3.dta"), 
                select = c(id,
                           Periode, 
                           Q26_1, Q26_2, Q26_4, Q26_5, 
                           q27_8, q27_9, Q29_1, 
                           Q30b_1, Q30b_2,
                           Q6_10, Q6_11,
                           Q34, 
                           kjonn, 
                           alder,
                           utdanning, 
                           persinnt12,
                           husinnt, 
                           vekt))



# __ Preparations ---------------------------------------------------------

# ___ Benefit reception: Self / Family ------------------------------------

# Self:
# NOTE: If a person received Unemployment benefits (Q26_1), 
#       Disability benefits (Q26_2), Social assistance (Q26_4), or 
#       SickPay/Sickness allowance (Q26_5), they are coded "benefit == 1"

SUPPA <- mutate(SUPPA, benefit = if_else(Q26_1 == 1 | Q26_2 == 1 | Q26_4 == 1 | Q26_5 == 1, 1, 0))
SUPPA$benefit <- if_else(is.na(SUPPA$Q26_1) | is.na(SUPPA$Q26_2) | is.na(SUPPA$Q26_4) | is.na(SUPPA$Q26_5), NA, SUPPA$benefit)
#table(SUPPA$benefit, exclude = NULL)



# ___ Lagging the income variable -----------------------------------------
SUPPA <- SUPPA %>%                           
  group_by(id) %>%
  dplyr::mutate(incomelag = dplyr::lag(persinnt12, n = 1, default = NA))



# ___ Income difference: difference: income - income_lagged ---------------
SUPPA <- mutate(SUPPA, incomediff = persinnt12 - incomelag)
#hist(SUPPA$incomediff)
#table(SUPPA$incomediff)



# ____ Dummies: Income increase and decrease ------------------------------
SUPPA <- dplyr::mutate(SUPPA, inc_increase = if_else(incomediff > 0, 1, 0))
SUPPA <- dplyr::mutate(SUPPA, inc_decrease = if_else(incomediff < 0, 1, 0))



# ____ Income groups ------------------------------------------------------
# Creating chronological sequence of income changes
income.seq <- 
  SUPPA %>%
  dplyr::select(id, Periode, persinnt12) %>%
  pivot_wider(names_from = Periode, values_from = persinnt12) %>%
  rename(income_1 = `1`,
         income_2 = `2`,
         income_3 = `3`)


# Group coding
SUPPA <-
  income.seq %>% 
  dplyr::mutate(income_grp = if_else(income_1 <= 3, "Low income",
                                     if_else(income_1 >=7, "High income", "Medium income"))) %>%
  dplyr::mutate(income_grp = if_else(!is.na(income_grp), income_grp,
                                     if_else(income_2 <= 3, "Low income",
                                             if_else(income_2 >=7, "High income", "Medium income")))) %>%
  dplyr::mutate(income_grp = if_else(!is.na(income_grp), income_grp,
                                     if_else(income_3 <= 3, "Low income",
                                             if_else(income_3 >=7, "High income", "Medium income")))) %>%
  dplyr::select(id, income_grp) %>%
  unique() %>%
  full_join(SUPPA, by = "id")

# Cleanup
rm(income.seq)



# ___ Female --------------------------------------------------------------
SUPPA <- mutate(SUPPA, male = if_else(kjonn == 2, 0, 1))



# ___ Removing NAs --------------------------------------------------------
SUPPA$Q30b_1  <- na_if(SUPPA$Q30b_1, 5)
SUPPA$Q30b_2  <- na_if(SUPPA$Q30b_2, 5)
SUPPA$Q34     <- na_if(SUPPA$Q34, 12)
SUPPA$husinnt <- na_if(SUPPA$husinnt, 9)
SUPPA$q27_8   <- na_if(SUPPA$q27_8, 8)
SUPPA$q27_9   <- na_if(SUPPA$q27_9, 8)
SUPPA$Q29_1   <- na_if(SUPPA$Q29_1, 12)
SUPPA$Q6_10   <- na_if(SUPPA$Q6_10, 5)
SUPPA$Q6_11   <- na_if(SUPPA$Q6_11, 5)




# ANALYSIS ----------------------------------------------------------------

# _ FIGURE 1 (Barplot): Answer distribution of responses in DVs -----------

# __ Generating data ------------------------------------------------------
barplot.pooled <-
  rbind(
    # Benefit abuse: Pooled all waves
    data.frame(
      SUPPA %>%
        ungroup %>%
        count(Answer = Q30b_1) %>%
        drop_na() %>%
        mutate(Group_total = sum(n)) %>%
        mutate(Share = round((n/Group_total)*100, 1)) %>%
        mutate(Type = c("Pooled"), .before = Answer) %>%
        mutate(DV = c("Benefit abuse"), .before = Answer)
    ),
    
    # Tax evasion: Pooled all waves
    data.frame(
      SUPPA %>%
        ungroup %>%
        count(Answer = Q30b_2) %>%
        drop_na() %>%
        mutate(Group_total = sum(n)) %>%
        mutate(Share = round((n/Group_total)*100, 1)) %>%
        mutate(Type = c("Pooled"), .before = Answer) %>%
        mutate(DV = c("Tax evasion"), .before = Answer)
    )
  )




# __ Figure 1 -------------------------------------------------------------
barplot.pool.fig <-
  annotate_figure(
    ggarrange(
      # Benefit abuse: Pooled
      ggplot(filter(barplot.pooled, DV == "Benefit abuse"), 
             aes(x = Answer, y = Share)) + 
        
        # Grid
        geom_hline(yintercept = seq(5, 45, 5), color = "grey90") +
        
        # Bars
        geom_bar(stat = "identity",
                 position = position_dodge(),
                 color = "black",
                 alpha = 0.66) +
        geom_text(aes(label = paste0(Share, "%\nN=", n)), 
                  color = "black", 
                  size = 9/.pt, 
                  position = position_dodge(0.9), 
                  alpha = 0.85, 
                  vjust = -0.1, 
                  lineheight = 0.8) +
        
        # Header
        geom_rect(xmin = 0.35, xmax = 4.65, ymin = 50, ymax = 62.5, colour = "black", fill = "grey90") +
        annotate("text", x = 2.5, y = 57.5, hjust = 0.5, vjust = -0.15, label = "Benefit abuse", size = 12/.pt, fontface = "bold") +
        annotate("text", x = 2.5, y = 57.5, hjust = 0.5, vjust =  1.25, label = "To remain at home and receive sick allowance\neven though they could have been working", size = 8.5/.pt, lineheight = 0.8) +
        
        # Specs
        scale_y_continuous(breaks = seq(0, 50, 10),
                           labels = c("0%", "10%", "20%", "30%", "40%", "50%")) +
        scale_x_continuous(breaks = c(1, 2, 3, 4), 
                           labels = c("Very\nuncommon", "Rather\nuncommon", "Rather\ncommon", "Very\ncommon")) +
        labs(x = NULL, y = NULL, fill = NULL) +
        coord_cartesian(xlim = c(0.35, 4.65), ylim = c(0, 62.5), expand = FALSE, clip = "off") +
        theme(text = element_text(color = "black"),
              legend.spacing.y = unit(0.25, "cm"),
              legend.key.size = unit(0.65, "cm"),
              panel.grid = element_blank(), 
              panel.background = element_rect(fill = "white"),
              plot.background = element_rect(fill = "white"),
              panel.border = element_rect(fill = "transparent"),
              axis.ticks = element_blank(),
              plot.margin = unit(c(0.2, 0.02, 0.2, 0.2),"cm")),
      
      
      # Tax evasion: Pooled
      ggplot(filter(barplot.pooled, DV == "Tax evasion"), 
             aes(x = Answer, y = Share)) + 
        
        # Grid
        geom_hline(yintercept = seq(5, 45, 5), color = "grey90") +
        
        # Bars
        geom_bar(stat = "identity",
                 position = position_dodge(),
                 color = "black",
                 alpha = 0.66) +
        geom_text(aes(label = paste0(Share, "%\nN=", n)), 
                  color = "black", 
                  size = 9/.pt, 
                  position = position_dodge(0.9), 
                  alpha = 0.85, 
                  vjust = -0.1, 
                  lineheight = 0.8) +
        
        # Header
        geom_rect(xmin = 0.35, xmax = 4.65, ymin = 50, ymax = 62.5, colour = "black", fill = "grey90") +
        annotate("text", x = 2.5, y = 57.5, hjust = 0.5, vjust = -0.15, label = "Tax evasion", size = 12/.pt, fontface = "bold") +
        annotate("text", x = 2.5, y = 57.5, hjust = 0.5, vjust =  1.25, label = "To misreport information about their economic\nsituation in order to avoid taxes or fees", size = 8.5/.pt, lineheight = 0.8) +
        
        # Specs
        scale_y_continuous(breaks = seq(0, 50, 10),
                           labels = c("0%", "10%", "20%", "30%", "40%", "50%")) +
        scale_x_continuous(breaks = c(1, 2, 3, 4), 
                           labels = c("Very\nuncommon", "Rather\nuncommon", "Rather\ncommon", "Very\ncommon")) +
        labs(x = NULL, y = NULL, fill = NULL) +
        coord_cartesian(xlim = c(0.35, 4.65), ylim = c(0, 62.5), expand = FALSE, clip = "off") +
        theme(text = element_text(color = "black"),
              axis.text.y = element_blank(),
              legend.spacing.y = unit(0.25, "cm"),
              legend.key.size = unit(0.65, "cm"),
              panel.grid = element_blank(), 
              panel.background = element_rect(fill = "white"),
              plot.background = element_rect(fill = "white"),
              panel.border = element_rect(fill = "transparent"),
              axis.ticks = element_blank(),
              plot.margin = unit(c(0.2, 0.85, 0.2, 0.2),"cm")),
      
      ncol = 2), 
    
    top = text_grob("How common or uncommon would you say it is\nfor people living in Norway to do the following?\n", size = 12, lineheight = 0.8)) +
  theme(panel.background = element_rect(fill = "white"),
        plot.background = element_rect(fill = "white"),
        panel.border = element_rect(fill = "transparent"))


# Saving the output
ggsave(barplot.pool.fig, filename = "04-figures and models/Fig 1 - barplots per waves.png",
       width = 17.5, height = 10, units = "cm", dpi = 600)


# Interim cleanup
rm(barplot.pooled, barplot.pool.fig)





# _ CROSS-SECTION: Between-group differences ------------------------------

# __ Table 1 --------------------------------------------------------------

# Models 1 & 2: Benefit abuse
ols.benefit <- feols(Q30b_1 ~ benefit +
                       sw0(
                         as.factor(Periode) + 
                           male + alder + Q34 + as.factor(persinnt12) + as.factor(utdanning)
                       ),
                     cluster = ~id, 
                     weights = ~vekt, 
                     data = na.omit(SUPPA[c("Q30b_1",
                                            "benefit", 
                                            "Periode", 
                                            "male", 
                                            "alder", 
                                            "Q34", 
                                            "persinnt12", 
                                            "utdanning", 
                                            "id", 
                                            "vekt")]))


# Models 3 & 4: Tax evasion
ols.tax <- feols(Q30b_2 ~ inc_increase + inc_decrease +
                   sw0(
                     as.factor(Periode) + 
                       male + alder + Q34 + as.factor(persinnt12) + as.factor(utdanning)
                   ),
                 cluster = ~id, 
                 weights = ~vekt, 
                 data = na.omit(SUPPA[c("Q30b_2", 
                                        "inc_increase", 
                                        "inc_decrease", 
                                        "Periode", 
                                        "male", 
                                        "alder", 
                                        "Q34", 
                                        "persinnt12", 
                                        "utdanning", 
                                        "id", 
                                        "vekt")]))


# Summary
etable(ols.benefit, ols.tax)



# Exporting the table
write_tableHTML(
  tableHTML(etable(ols.benefit, ols.tax, tex = FALSE, digits = 3)), 
  file = "04-figures and models/Tab 1 - between-group differences.html")



# __ Figure 2 -------------------------------------------------------------

# Extracting the marginal effects
margins <-
  rbind(
    data.frame(topic = "Benefit abuse",
               group = "Neither I, nor a\nfamily member\nreceived",
               estimate = 0,
               conf.low = NA,
               conf.high = NA,
               conf.low.1 = NA,
               conf.high.1 = NA),
    data.frame(topic = c("Benefit abuse"),
               group = c("Received myself"),
               broom::tidy(ols.benefit[[2]], conf.int = TRUE, conf.level = .95)[2, 2],
               broom::tidy(ols.benefit[[2]], conf.int = TRUE, conf.level = .95)[2, 6],
               broom::tidy(ols.benefit[[2]], conf.int = TRUE, conf.level = .95)[2, 7],
               broom::tidy(ols.benefit[[2]], conf.int = TRUE, conf.level = .90)[2, 6],
               broom::tidy(ols.benefit[[2]], conf.int = TRUE, conf.level = .90)[2, 7]),
    data.frame(topic = "Tax evasion",
               group = "No change",
               estimate = 0,
               conf.low = NA,
               conf.high = NA,
               conf.low.1 = NA,
               conf.high.1 = NA),
    data.frame(topic = c("Tax evasion"),
               group = c("Income\nincreased",
                         "Income\ndecreased"),
               broom::tidy(ols.tax[[2]], conf.int = TRUE, conf.level = .95)[2:3, 2],
               broom::tidy(ols.tax[[2]], conf.int = TRUE, conf.level = .95)[2:3, 6],
               broom::tidy(ols.tax[[2]], conf.int = TRUE, conf.level = .95)[2:3, 7],
               broom::tidy(ols.tax[[2]], conf.int = TRUE, conf.level = .90)[2:3, 6],
               broom::tidy(ols.tax[[2]], conf.int = TRUE, conf.level = .90)[2:3, 7])
  )



# Figure
ggpubr::ggarrange(
  
  # Benefit abuse
  ggplot(subset(margins, topic == "Benefit abuse"), aes(x = estimate, y = group)) +
    geom_vline(xintercept = 0, colour = "grey60", linewidth = 1, alpha = 0.5) +
    geom_errorbar(aes(xmin = conf.low,   xmax = conf.high), width = 0, linewidth = 0.8, alpha = 0.6) +
    geom_errorbar(aes(xmin = conf.low.1, xmax = conf.high.1), width = 0, linewidth = 1.25) +
    geom_point(size = 2) +
    labs(y = NULL, x = NULL, title = NULL) +
    scale_y_discrete(limits = rev(c("Received myself",
                                    "Neither I, nor a\nfamily member\nreceived")),
                     labels = rev(c("Received",
                                    "Not received"))) +
    geom_rect(ymin = 2.5, ymax = 3, xmin = -0.5, xmax = 0.5, colour = "black", fill = "grey90") +
    annotate("text", x = 0, y = 2.75, hjust = 0.5, vjust = 0.5, label = "Benefit abuse", size = 11/.pt) +
    coord_cartesian(xlim = c(-0.5, 0.5), ylim = c(0.5, 3), expand = FALSE, clip = "off") +
    scale_x_continuous(breaks = seq(-0.5, 0.5, 0.1)) +
    theme(panel.grid.major.y = element_blank(), 
          panel.grid.major.x = element_line(colour = "gray90"), 
          panel.grid.minor.x = element_blank(), 
          panel.background = element_rect(fill = "white"),
          plot.background = element_rect(fill = "white"),
          panel.border = element_rect(fill = "transparent"),
          axis.ticks = element_blank(),
          axis.text.x = element_blank(),
          plot.margin = unit(c(0.5, 0.5, 0.15, 0.5),"cm")),
  
  # Tax evasion
  ggplot(subset(margins, topic == "Tax evasion"), aes(x = estimate, y = group)) +
    geom_vline(xintercept = 0, colour = "grey60", linewidth = 1, alpha = 0.5) +
    geom_errorbar(aes(xmin = conf.low,   xmax = conf.high), width = 0, linewidth = 0.8, alpha = 0.6) +
    geom_errorbar(aes(xmin = conf.low.1, xmax = conf.high.1), width = 0, linewidth = 1.25) +
    geom_point(size = 2) +
    labs(y = NULL, x = "Perception of system abuse: Group differences", title = NULL) +
    scale_y_discrete(limits = rev(c("Income\nincreased",
                                    "Income\ndecreased", 
                                    "No change"))) +
    geom_rect(ymin = 3.5, ymax = 4, xmin = -0.5, xmax = 0.5, colour = "black", fill = "grey90") +
    annotate("text", x = 0, y = 3.75, hjust = 0.5, vjust = 0.5, label = "Tax evasion", size = 11/.pt) +
    coord_cartesian(xlim = c(-0.5, 0.5), ylim = c(0.5, 4), expand = FALSE, clip = "off") +
    scale_x_continuous(breaks = seq(-0.5, 0.5, 0.1)) +
    theme(panel.grid.major.y = element_blank(), 
          panel.grid.major.x = element_line(colour = "gray90"), 
          panel.grid.minor.x = element_blank(), 
          panel.background = element_rect(fill = "white"),
          plot.background = element_rect(fill = "white"),
          panel.border = element_rect(fill = "transparent"),
          axis.ticks = element_blank(),
          plot.margin = unit(c(0.15, 0.5, 0.5, 0.5),"cm")),
  
  nrow = 2, heights = c(7, 11),
  align = "v")


# Saving the output
ggsave(filename = "04-figures and models/Fig 2 - ols - group comparison.png", width = 15, height = 10, units = "cm", dpi = 600)


# Cleanup
rm(ols.benefit, ols.tax, margins)





# _ FIXED EFFECTS: Within-individual effects ------------------------------

# __ Table 2 --------------------------------------------------------------

# ___ SICKNESS BENEFITS: Models -------------------------------------------

# BASE model for sickness benefits (i.e., benefit abuse)
fix.benefit.base <- feols(Q30b_1 ~ benefit | id + Periode,  
                          cluster = ~id, 
                          weights = ~vekt, 
                          data = SUPPA)

# EFFECT HETEROGENEITY for sickness benefits (based on initial income)
fix.benefit.inc_grp <- feols(Q30b_1 ~ benefit * income_grp | id + Periode,  
                            cluster = ~id, 
                            weights = ~vekt, 
                            data = SUPPA)





# ___ INCOME CHANGE: Models -----------------------------------------------

# BASE model for income change (i.e., tax evasion)
fix.tax.base <- feols(Q30b_2 ~ persinnt12 | id + Periode,  
                      cluster = ~id, 
                      weights = ~vekt, 
                      data = SUPPA)

# EFFECT HETEROGENEITY for income change (based on initial income)
fix.tax.inc_grp <- feols(Q30b_2 ~ persinnt12 * income_grp | id + Periode,  
                          cluster = ~id, 
                          weights = ~vekt, 
                          data = SUPPA)



# Output
etable(fix.benefit.base, fix.benefit.inc_grp, fix.tax.base, fix.tax.inc_grp)




# Exporting the table
write_tableHTML(
  tableHTML(etable(fix.benefit.base, fix.benefit.inc_grp, fix.tax.base,fix.tax.inc_grp, tex = FALSE, digits = 3)), 
  file = "04-figures and models/Tab 2 - TWFE models.html")





# __ Figure 3 -------------------------------------------------------------

# Re-estimating the models with interactions to get straightforward coefficients
fix.benefit.inc_grp <- feols(Q30b_1 ~ benefit : income_grp | id + Periode,  
                             cluster = ~id, 
                             weights = ~vekt, 
                             data = SUPPA)

fix.tax.inc_grp <- feols(Q30b_2 ~ persinnt12 : income_grp | id + Periode,  
                         cluster = ~id, 
                         weights = ~vekt, 
                         data = SUPPA)




# Extracting the marginal effects
margins <-
  rbind(
    data.frame(topic = "Benefit abuse",
               group = c("General",
                         "High income",
                         "Low income",
                         "Medium income"),
               time = c("Time 0"),
               estimate = 0,
               conf.low = NA,
               conf.high = NA),
    data.frame(topic = c("Benefit abuse"),
               group = "General",
               time = c("Time 1"),
               broom::tidy(fix.benefit.base, conf.int = TRUE, conf.level = .95)[1, c(2, 6:7)]),
    data.frame(topic = c("Benefit abuse"),
               group = c("High income",
                         "Low income",
                         "Medium income"),
               time = c("Time 1"),
               broom::tidy(fix.benefit.inc_grp, conf.int = TRUE, conf.level = .95)[1:3, c(2, 6:7)]),
    data.frame(topic = "Tax evasion",
               group = c("General",
                         "Low income",
                         "Medium income",
                         "High income"),
               time = c("Time 0"),
               estimate = 0,
               conf.low = NA,
               conf.high = NA),
    data.frame(topic = c("Tax evasion"),
               group = "General",
               time = c("Time 1"),
               broom::tidy(fix.tax.base, conf.int = TRUE, conf.level = .95)[1, c(2, 6:7)]),
    data.frame(topic = c("Tax evasion"),
               group = c("High income",
                         "Low income",
                         "Medium income"),
               time = c("Time 1"),
               broom::tidy(fix.tax.inc_grp, conf.int = TRUE, conf.level = .95)[1:3, c(2, 6:7)])
  )



# Figure
annotate_figure(
  cowplot::plot_grid(
    
    # Benefit abuse: General
    ggplot(filter(margins, topic == "Benefit abuse" & group == "General"), 
           aes(x = time, y = estimate, group = group)) +
      geom_hline(yintercept = c(seq(-0.3, -0.05, 0.05), seq(0.05, 0.3, 0.05)), colour = "grey90", linewidth = 0.75, alpha = 0.5) +
      geom_hline(yintercept = 0, colour = "grey60", linewidth = 1, alpha = 0.5) +
      geom_point() +
      geom_line() +
      geom_errorbar(aes(ymin = conf.low, ymax = conf.high), width = 0.025) +
      labs(y = NULL, x = NULL, title = NULL) +
      scale_y_continuous(breaks = seq(-0.3, 0.3, 0.1),
                         labels = c("-0.3", "-0.2", "-0.1", "0.0", "0.1", "0.2", "0.3")) +
      geom_rect(ymin = 0.3, ymax = 0.4, xmin = 0.5, xmax = 2.5, colour = "black", fill = "grey90") +
      annotate("text", x = 1.5, y = 0.35, hjust = 0.5, vjust = 0.5, label = "Benefit abuse", size = 11/.pt) +
      coord_cartesian(xlim = c(0.5, 2.5), ylim = c(-0.3, 0.4), expand = FALSE, clip = "off") +
      theme(panel.grid = element_blank(), 
            panel.background = element_rect(fill = "white"),
            plot.background = element_rect(fill = "white"),
            panel.border = element_rect(fill = "transparent"),
            axis.ticks = element_blank(),
            axis.text.x = element_blank(),
            plot.margin = unit(c(0.5, 0.5, 0.5, 0.5),"cm")),
    
    # Tax evasion: General
    ggplot(filter(margins, topic == "Tax evasion" & group == "General"), 
           aes(x = time, y = estimate, group = group)) +
      geom_hline(yintercept = c(seq(-0.3, -0.05, 0.05), seq(0.05, 0.3, 0.05)), colour = "grey90", linewidth = 0.75, alpha = 0.5) +
      geom_hline(yintercept = 0, colour = "grey60", linewidth = 1, alpha = 0.5) +
      geom_point() +
      geom_line() +
      geom_errorbar(aes(ymin = conf.low, ymax = conf.high), width = 0.025) +
      labs(y = NULL, x = NULL, title = NULL) +
      scale_y_continuous(breaks = seq(-0.3, 0.3, 0.1),
                         labels = c("-0.3", "-0.2", "-0.1", "0.0", "0.1", "0.2", "0.3")) +
      geom_rect(ymin = 0.3, ymax = 0.4, xmin = 0.5, xmax = 2.5, colour = "black", fill = "grey90") +
      annotate("text", x = 1.5, y = 0.35, hjust = 0.5, vjust = 0.5, label = "Tax evasion", size = 11/.pt) +
      coord_cartesian(xlim = c(0.5, 2.5), ylim = c(-0.3, 0.4), expand = FALSE, clip = "off") +
      theme(panel.grid = element_blank(), 
            panel.background = element_rect(fill = "white"),
            plot.background = element_rect(fill = "white"),
            panel.border = element_rect(fill = "transparent"),
            axis.ticks = element_blank(),
            axis.text = element_blank(),
            plot.margin = unit(c(0.5, 0.5, 0.5, 0.5),"cm")),
    
    NULL,
    
    # Benefit abuse: Heterogeneity
    ggplot(filter(margins, topic == "Benefit abuse" & group != "General"), 
           aes(x = time, y = estimate, group = group, color = group, shape = group)) +
      geom_hline(yintercept = c(seq(-0.3, -0.05, 0.05), seq(0.05, 0.3, 0.05)), colour = "grey90", linewidth = 0.75, alpha = 0.5) +
      geom_hline(yintercept = 0, colour = "grey60", linewidth = 1, alpha = 0.5) +
      geom_point(position = position_dodge(width = 0.3)) +
      geom_line(position = position_dodge(width = 0.3)) +
      geom_errorbar(aes(ymin = conf.low, ymax = conf.high), width = 0.025, position = position_dodge(width = 0.3)) +
      labs(y = NULL, x = NULL, title = NULL, color = "Initial income level", shape = "Initial income level") +
      scale_y_continuous(breaks = seq(-0.3, 0.3, 0.1),
                         labels = c("-0.3", "-0.2", "-0.1", "0.0", "0.1", "0.2", "0.3")) +
      scale_color_manual(breaks = c("High income", "Medium income", "Low income"),
                         labels = c("High [7-9]", "Medium [4-6]", "Low [1-3]"),
                         values = c(brewer.pal(n = 9, "Greys")[c(4, 6, 9)])) +
      scale_shape_manual(breaks = c("High income", "Medium income", "Low income"),
                         labels = c("High [7-9]", "Medium [4-6]", "Low [1-3]"),
                         values = c(16, 17, 15)) +
      scale_x_discrete(breaks = c("Time 0", "Time 1"),
                       labels = c("<span style='font-size:9pt'>Time <sub>*t=0*</sub></span><br/>",
                                  "<span style='font-size:9pt'>Time <sub>*t=+1*</sub></span><br/>")) +
      coord_cartesian(xlim = c(0.5, 2.5), ylim = c(-0.3, 0.3), expand = FALSE, clip = "off") +
      theme(panel.grid = element_blank(), 
            panel.background = element_rect(fill = "white"),
            plot.background = element_rect(fill = "white"),
            panel.border = element_rect(fill = "transparent"),
            legend.spacing.y = unit(0.2, 'cm'),
            legend.key = element_rect(linewidth = 10, fill = "white"),
            legend.key.size = unit(0.5, 'cm'),
            legend.title = element_text(face = "bold"),
            legend.position = "none",
            axis.ticks = element_blank(),
            axis.text.x = element_markdown(hjust = 0.5),
            plot.margin = unit(c(0.5, 0.5, 0.5, 0.5),"cm")),
    
    # Tax evasion: Heterogeneity
    ggplot(filter(margins, topic == "Tax evasion" & group != "General"), 
           aes(x = time, y = estimate, group = group, color = group, shape = group)) +
      geom_hline(yintercept = c(seq(-0.3, -0.05, 0.05), seq(0.05, 0.3, 0.05)), colour = "grey90", linewidth = 0.75, alpha = 0.5) +
      geom_hline(yintercept = 0, colour = "grey60", linewidth = 1, alpha = 0.5) +
      geom_point(position = position_dodge(width = 0.3)) +
      geom_line(position = position_dodge(width = 0.3)) +
      geom_errorbar(aes(ymin = conf.low, ymax = conf.high), width = 0.025, position = position_dodge(width = 0.3)) +
      labs(y = NULL, x = NULL, title = NULL, color = "Initial income level", shape = "Initial income level") +
      scale_y_continuous(breaks = seq(-0.3, 0.3, 0.1),
                         labels = c("-0.3", "-0.2", "-0.1", "0.0", "0.1", "0.2", "0.3")) +
      scale_color_manual(breaks = c("High income", "Medium income", "Low income"),
                         labels = c("High [7-9]", "Medium [4-6]", "Low [1-3]"),
                         values = c(brewer.pal(n = 9, "Greys")[c(4, 6, 9)])) +
      scale_shape_manual(breaks = c("High income", "Medium income", "Low income"),
                         labels = c("High [7-9]", "Medium [4-6]", "Low [1-3]"),
                         values = c(16, 17, 15)) +
      scale_x_discrete(breaks = c("Time 0", "Time 1"),
                       labels = c("<span style='font-size:9pt'>Time <sub>*t=0*</sub></span><br/>",
                                  "<span style='font-size:9pt'>Time <sub>*t=+1*</sub></span><br/>")) +
      coord_cartesian(xlim = c(0.5, 2.5), ylim = c(-0.3, 0.3), expand = FALSE, clip = "off") +
      theme(panel.grid = element_blank(), 
            panel.background = element_rect(fill = "white"),
            plot.background = element_rect(fill = "white"),
            panel.border = element_rect(fill = "transparent"),
            legend.spacing.y = unit(0.2, 'cm'),
            legend.key = element_rect(linewidth = 10, fill = "white"),
            legend.key.size = unit(0.5, 'cm'),
            legend.title = element_text(face = "bold"),
            legend.position = "none",
            axis.ticks = element_blank(),
            axis.text.y = element_blank(),
            axis.text.x = element_markdown(hjust = 0.5),
            plot.margin = unit(c(0.5, 0.5, 0.5, 0.5),"cm")),
    
    get_legend(
      ggplot(filter(margins, topic == "Benefit abuse" & group != "General"), 
             aes(x = time, y = estimate, group = group, color = group, shape = group)) +
        geom_hline(yintercept = c(seq(-0.3, -0.05, 0.05), seq(0.05, 0.3, 0.05)), colour = "grey90", linewidth = 0.75, alpha = 0.5) +
        geom_hline(yintercept = 0, colour = "grey60", linewidth = 1.25, alpha = 0.5) +
        geom_point(position = position_dodge(width = 0.3)) +
        geom_line(position = position_dodge(width = 0.3)) +
        geom_errorbar(aes(ymin = conf.low, ymax = conf.high), width = 0.025, position = position_dodge(width = 0.3)) +
        labs(y = NULL, x = NULL, title = NULL, color = "Initial income level", shape = "Initial income level") +
        scale_y_continuous(breaks = seq(-0.3, 0.3, 0.1),
                           labels = c("-0.3", "-0.2", "-0.1", "0.0", "0.1", "0.2", "0.3")) +
        scale_color_manual(breaks = c("High income", "Medium income", "Low income"),
                           labels = c("<span style='font-size:9pt'>High</span><br/><span style='font-size:7pt'>Above 700,000 NOK</span><br/>", 
                                      "<span style='font-size:9pt'>Medium</span><br/><span style='font-size:7pt'>400,000 to 699,999 NOK</span><br/>", 
                                      "<span style='font-size:9pt'>Low</span><br/><span style='font-size:7pt'>Below 399,999 NOK</span>"),
                           values = c(brewer.pal(n = 9, "Greys")[c(4, 6, 9)])) +
        scale_shape_manual(breaks = c("High income", "Medium income", "Low income"),
                           labels = c("<span style='font-size:9pt'>High</span><br/><span style='font-size:7pt'>Above 700,000 NOK</span><br/>", 
                                      "<span style='font-size:9pt'>Medium</span><br/><span style='font-size:7pt'>400,000 to 699,999 NOK</span><br/>", 
                                      "<span style='font-size:9pt'>Low</span><br/><span style='font-size:7pt'>Below 399,999 NOK</span>"),
                           values = c(16, 17, 15)) +
        scale_x_discrete(breaks = c("Time 0", "Time 1"),
                         labels = c("<span style='font-size:9pt'>Time <sub>*t=0*</sub></span><br/>",
                                    "<span style='font-size:9pt'>Time <sub>*t=+1*</sub></span><br/>")) +
        coord_cartesian(xlim = c(0.5, 2.5), ylim = c(-0.3, 0.3), expand = FALSE, clip = "off") +
        theme(panel.grid = element_blank(), 
              panel.background = element_rect(fill = "white"),
              plot.background = element_rect(fill = "white"),
              panel.border = element_rect(fill = "transparent"),
              legend.key = element_rect(linewidth = 12, fill = "white"),
              legend.key.size = unit(0.75, 'cm'),
              legend.title = element_text(face = "bold", size = 26/.pt),
              legend.text = element_markdown(),
              legend.position = c(0.5, 0.57),
              legend.background = element_rect(fill = "transparent"),
              axis.ticks = element_blank(),
              axis.text.x = element_markdown(hjust = 0.5),
              plot.margin = unit(c(0.5, 0.5, 0.5, 0.75),"cm"))),
    
    ncol = 3, rel_widths = c(10, 9, 4),
    nrow = 2,
    align = "none",
    common.legend = TRUE, legend = "right"),
  
  left =  text_grob("Difference in the perception of system abuse", rot = 90, face = "bold")
) + theme(panel.background = element_rect(fill = "white", color = "white"),
          plot.background = element_rect(fill = "white", color = "white"),
          panel.border = element_rect(fill = "transparent", color = "white"),
          plot.margin = unit(c(0, 0.75, 0, 0.15),"cm"))


# Saving the output
ggsave(filename = "04-figures and models/Fig 3 - fixed effects.png", width = 17.5, height = 15, units = "cm", dpi = 600)



# Cleanup
rm(fix.benefit, fix.tax, margins)
rm(fix.benefit.base, fix.benefit.inc_grp, fix.tax.base, fix.tax.inc_grp)





# ONLINE APPENDIX ---------------------------------------------------------

# _ DESCRIPTIVE STATISTICS ------------------------------------------------

# __ Table S3-S4: Descriptive statistics ----------------------------------

# Wave 1
dfSummary(
  subset(SUPPA, Periode == 1)[c("Q30b_1",
                                "Q30b_2",
                                "benefit", 
                                "inc_increase", 
                                "inc_decrease", 
                                "Periode", 
                                "male", 
                                "alder", 
                                "Q34", 
                                "persinnt12", 
                                "utdanning")],
  varnumbers = FALSE,
  labels.col = FALSE,
  graph.col = FALSE,
  display.labels = FALSE,
  na.col = FALSE)



# Wave 2
dfSummary(
  subset(SUPPA, Periode == 2)[c("Q30b_1",
                                "Q30b_2",
                                "benefit", 
                                "inc_increase", 
                                "inc_decrease", 
                                "Periode", 
                                "male", 
                                "alder", 
                                "Q34", 
                                "persinnt12", 
                                "utdanning")],
  varnumbers = FALSE,
  labels.col = FALSE,
  graph.col = FALSE,
  display.labels = FALSE,
  na.col = FALSE)



# Wave 3
dfSummary(
  subset(SUPPA, Periode == 3)[c("Q30b_1",
                                "Q30b_2",
                                "benefit", 
                                "inc_increase", 
                                "inc_decrease", 
                                "Periode", 
                                "male", 
                                "alder", 
                                "Q34", 
                                "persinnt12", 
                                "utdanning")],
  varnumbers = FALSE,
  labels.col = FALSE,
  graph.col = FALSE,
  display.labels = FALSE,
  na.col = FALSE)



# Total
dfSummary(
  SUPPA[c("Q30b_1",
          "Q30b_2",
          "benefit", 
          "inc_increase", 
          "inc_decrease", 
          "Periode", 
          "male", 
          "alder", 
          "Q34", 
          "persinnt12", 
          "utdanning")],
  varnumbers = FALSE,
  labels.col = FALSE,
  graph.col = FALSE,
  display.labels = FALSE,
  na.col = FALSE)




# _ FIGURE S1 (Barplots): Answer distribution in DVs per wave -------------

# _ Generating data -------------------------------------------------------

barplot.waves <-
  rbind(
    # Benefit abuse: Wave 1
    data.frame(
      SUPPA %>%
        ungroup %>%
        filter(Periode == 1) %>%
        count(Answer = Q30b_1) %>%
        drop_na() %>%
        mutate(Group_total = sum(n)) %>%
        mutate(Share = round((n/Group_total)*100, 1)) %>%
        mutate(Type = c("Ind. waves"), .before = Answer) %>%
        mutate(DV = c("Benefit abuse"), .before = Answer) %>%
        mutate(Wave = c("Wave 1"), .before = Answer)
    ),
    
    # Benefit abuse: Wave 2
    data.frame(
      SUPPA %>%
        ungroup %>%
        filter(Periode == 2) %>%
        count(Answer = Q30b_1) %>%
        drop_na() %>%
        mutate(Group_total = sum(n)) %>%
        mutate(Share = round((n/Group_total)*100, 1)) %>%
        mutate(Type = c("Ind. waves"), .before = Answer) %>%
        mutate(DV = c("Benefit abuse"), .before = Answer) %>%
        mutate(Wave = c("Wave 2"), .before = Answer)
    ),
    
    # Benefit abuse: Wave 3
    data.frame(
      SUPPA %>%
        ungroup %>%
        filter(Periode == 3) %>%
        count(Answer = Q30b_1) %>%
        drop_na() %>%
        mutate(Group_total = sum(n)) %>%
        mutate(Share = round((n/Group_total)*100, 1)) %>%
        mutate(Type = c("Ind. waves"), .before = Answer) %>%
        mutate(DV = c("Benefit abuse"), .before = Answer) %>%
        mutate(Wave = c("Wave 3"), .before = Answer)
    ),
    
    # Tax evasion: Wave 1
    data.frame(
      SUPPA %>%
        ungroup %>%
        filter(Periode == 1) %>%
        count(Answer = Q30b_2) %>%
        drop_na() %>%
        mutate(Group_total = sum(n)) %>%
        mutate(Share = round((n/Group_total)*100, 1)) %>%
        mutate(Type = c("Ind. waves"), .before = Answer) %>%
        mutate(DV = c("Tax evasion"), .before = Answer) %>%
        mutate(Wave = c("Wave 1"), .before = Answer)
    ),
    
    # Tax evasion: Wave 2
    data.frame(
      SUPPA %>%
        ungroup %>%
        filter(Periode == 2) %>%
        count(Answer = Q30b_2) %>%
        drop_na() %>%
        mutate(Group_total = sum(n)) %>%
        mutate(Share = round((n/Group_total)*100, 1)) %>%
        mutate(Type = c("Ind. waves"), .before = Answer) %>%
        mutate(DV = c("Tax evasion"), .before = Answer) %>%
        mutate(Wave = c("Wave 2"), .before = Answer)
    ),
    
    # Tax evasion: Wave 3
    data.frame(
      SUPPA %>%
        ungroup %>%
        filter(Periode == 3) %>%
        count(Answer = Q30b_2) %>%
        drop_na() %>%
        mutate(Group_total = sum(n)) %>%
        mutate(Share = round((n/Group_total)*100, 1)) %>%
        mutate(Type = c("Ind. waves"), .before = Answer) %>%
        mutate(DV = c("Tax evasion"), .before = Answer) %>%
        mutate(Wave = c("Wave 3"), .before = Answer)
    )
  )




# __ Figure S1 ------------------------------------------------------------

barplot.waves.fig <-
  ggarrange(
    
    # Benefit abuse: Wave 1
    ggplot(filter(barplot.waves, Wave == "Wave 1" & DV == "Benefit abuse"), 
           aes(x = Answer, y = Share)) + 
      # Grid
      geom_hline(yintercept = seq(5, 45, 5), color = "grey90") +
      
      # Bars
      geom_bar(stat = "identity",
               position = position_dodge(),
               color = "black",
               alpha = 0.66) +
      geom_text(aes(label = paste0(Share, "%\nN=", n)), 
                color = "black", 
                size = 9/.pt, 
                position = position_dodge(0.9), 
                alpha = 0.85, 
                vjust = -0.1, 
                lineheight = 0.8) +
      
      # Header
      geom_rect(xmin = 0.35, xmax = 4.65, ymin = 50, ymax = 60, colour = "black", fill = "grey90") +
      annotate("text", x = 2.5, y = 55, hjust = 0.5, vjust = 0.5, label = "Wave 1", size = 12/.pt) +
      
      # Specs
      scale_y_continuous(breaks = seq(0, 50, 10),
                         labels = c("0%", "10%", "20%", "30%", "40%", "50%")) +
      scale_x_continuous(breaks = c(1, 2, 3, 4), 
                         labels = c("Very\nuncommon", "Rather\nuncommon", "Rather\ncommon", "Very\ncommon")) +
      labs(x = NULL, y = NULL, fill = NULL) +
      coord_cartesian(xlim = c(0.35, 4.65), ylim = c(0, 60), expand = FALSE, clip = "off") +
      theme(text = element_text(color = "black"),
            axis.text.x = element_blank(),
            legend.spacing.y = unit(0.25, "cm"),
            legend.key.size = unit(0.65, "cm"),
            panel.grid = element_blank(), 
            panel.background = element_rect(fill = "white"),
            plot.background = element_rect(fill = "white"),
            panel.border = element_rect(fill = "transparent"),
            axis.ticks = element_blank(),
            plot.margin = unit(c(0.2, 0.02, 0.2, 0.2),"cm")),
    
    # Benefit abuse: Wave 2
    ggplot(filter(barplot.waves, Wave == "Wave 2" & DV == "Benefit abuse"), 
           aes(x = Answer, y = Share)) + 
      # Grid
      geom_hline(yintercept = seq(5, 45, 5), color = "grey90") +
      
      # Bars
      geom_bar(stat = "identity",
               position = position_dodge(),
               color = "black",
               alpha = 0.66) +
      geom_text(aes(label = paste0(Share, "%\nN=", n)), 
                color = "black", 
                size = 9/.pt, 
                position = position_dodge(0.9), 
                alpha = 0.85, 
                vjust = -0.1, 
                lineheight = 0.8) +
      
      # Header
      geom_rect(xmin = 0.35, xmax = 4.65, ymin = 50, ymax = 60, colour = "black", fill = "grey90") +
      annotate("text", x = 2.5, y = 55, hjust = 0.5, vjust = 0.5, label = "Wave 2", size = 12/.pt) +
      
      # Specs
      scale_y_continuous(breaks = seq(0, 50, 10),
                         labels = c("0%", "10%", "20%", "30%", "40%", "50%")) +
      scale_x_continuous(breaks = c(1, 2, 3, 4), 
                         labels = c("Very\nuncommon", "Rather\nuncommon", "Rather\ncommon", "Very\ncommon")) +
      labs(x = NULL, y = NULL, fill = NULL) +
      coord_cartesian(xlim = c(0.35, 4.65), ylim = c(0, 60), expand = FALSE, clip = "off") +
      theme(text = element_text(color = "black"),
            axis.text.x = element_blank(),
            axis.text.y = element_blank(),
            legend.spacing.y = unit(0.25, "cm"),
            legend.key.size = unit(0.65, "cm"),
            panel.grid = element_blank(), 
            panel.background = element_rect(fill = "white"),
            plot.background = element_rect(fill = "white"),
            panel.border = element_rect(fill = "transparent"),
            axis.ticks = element_blank(),
            plot.margin = unit(c(0.2, 0.02, 0.2, 0.2),"cm")),
    
    # Benefit abuse: Wave 3
    ggplot(filter(barplot.waves, Wave == "Wave 3" & DV == "Benefit abuse"), 
           aes(x = Answer, y = Share)) + 
      # Grid
      geom_hline(yintercept = seq(5, 45, 5), color = "grey90") +
      
      # Bars
      geom_bar(stat = "identity",
               position = position_dodge(),
               color = "black",
               alpha = 0.66) +
      geom_text(aes(label = paste0(Share, "%\nN=", n)), 
                color = "black", 
                size = 9/.pt, 
                position = position_dodge(0.9), 
                alpha = 0.85, 
                vjust = -0.1, 
                lineheight = 0.8) +
      
      # Header
      geom_rect(xmin = 0.35, xmax = 4.65, ymin = 50, ymax = 60, colour = "black", fill = "grey90") +
      annotate("text", x = 2.5, y = 55, hjust = 0.5, vjust = 0.5, label = "Wave 3", size = 12/.pt) +
      
      # Right panel
      geom_rect(xmin = 4.65, xmax = 5.15, ymin = 0, ymax = 50, colour = "black", fill = "grey90") +
      annotate("text", x = 4.9, y = 25, label = "Benefit abuse", angle = 270, size = 12/.pt) +
      
      # Specs
      scale_y_continuous(breaks = seq(0, 50, 10),
                         labels = c("0%", "10%", "20%", "30%", "40%", "50%")) +
      scale_x_continuous(breaks = c(1, 2, 3, 4), 
                         labels = c("Very\nuncommon", "Rather\nuncommon", "Rather\ncommon", "Very\ncommon")) +
      labs(x = NULL, y = NULL, fill = NULL) +
      coord_cartesian(xlim = c(0.35, 5.15), ylim = c(0, 60), expand = FALSE, clip = "off") +
      geom_hline(yintercept = 0) + geom_vline(xintercept = 0.35) +
      theme(text = element_text(color = "black"),
            axis.text.x = element_blank(),
            axis.text.y = element_blank(),
            legend.spacing.y = unit(0.25, "cm"),
            legend.key.size = unit(0.65, "cm"),
            panel.grid = element_blank(), 
            panel.background = element_rect(fill = "white"),
            plot.background = element_rect(fill = "white"),
            panel.border = element_rect(color = "transparent", fill = "transparent"),
            axis.ticks = element_blank(),
            plot.margin = unit(c(0.2, 0.02, 0.2, 0.2),"cm")),
    
    
    # Tax evasion: Wave 1
    ggplot(filter(barplot.waves, Wave == "Wave 1" & DV == "Tax evasion"), 
           aes(x = Answer, y = Share)) + 
      # Grid
      geom_hline(yintercept = seq(5, 45, 5), color = "grey90") +
      
      # Bars
      geom_bar(stat = "identity",
               position = position_dodge(),
               color = "black",
               alpha = 0.66) +
      geom_text(aes(label = paste0(Share, "%\nN=", n)), 
                color = "black", 
                size = 9/.pt, 
                position = position_dodge(0.9), 
                alpha = 0.85, 
                vjust = -0.1, 
                lineheight = 0.8) +
      
      # Specs
      scale_y_continuous(breaks = seq(0, 50, 10),
                         labels = c("0%", "10%", "20%", "30%", "40%", "50%")) +
      scale_x_continuous(breaks = c(1, 2, 3, 4), 
                         labels = c("Very\nuncommon", "Rather\nuncommon", "Rather\ncommon", "Very\ncommon")) +
      labs(x = NULL, y = NULL, fill = NULL) +
      coord_cartesian(xlim = c(0.35, 4.65), ylim = c(0, 50), expand = FALSE, clip = "off") +
      theme(text = element_text(color = "black"),
            legend.spacing.y = unit(0.25, "cm"),
            legend.key.size = unit(0.65, "cm"),
            panel.grid = element_blank(), 
            panel.background = element_rect(fill = "white"),
            plot.background = element_rect(fill = "white"),
            panel.border = element_rect(fill = "transparent"),
            axis.ticks = element_blank(),
            plot.margin = unit(c(0.2, 0.02, 0.2, 0.2),"cm")),
    
    # Tax evasion: Wave 2
    ggplot(filter(barplot.waves, Wave == "Wave 2" & DV == "Tax evasion"), 
           aes(x = Answer, y = Share)) + 
      # Grid
      geom_hline(yintercept = seq(5, 45, 5), color = "grey90") +
      
      # Bars
      geom_bar(stat = "identity",
               position = position_dodge(),
               color = "black",
               alpha = 0.66) +
      geom_text(aes(label = paste0(Share, "%\nN=", n)), 
                color = "black", 
                size = 9/.pt, 
                position = position_dodge(0.9), 
                alpha = 0.85, 
                vjust = -0.1, 
                lineheight = 0.8) +
      
      # Specs
      scale_y_continuous(breaks = seq(0, 50, 10),
                         labels = c("0%", "10%", "20%", "30%", "40%", "50%")) +
      scale_x_continuous(breaks = c(1, 2, 3, 4), 
                         labels = c("Very\nuncommon", "Rather\nuncommon", "Rather\ncommon", "Very\ncommon")) +
      labs(x = NULL, y = NULL, fill = NULL) +
      coord_cartesian(xlim = c(0.35, 4.65), ylim = c(0, 50), expand = FALSE, clip = "off") +
      theme(text = element_text(color = "black"),
            axis.text.y = element_blank(),
            legend.spacing.y = unit(0.25, "cm"),
            legend.key.size = unit(0.65, "cm"),
            panel.grid = element_blank(), 
            panel.background = element_rect(fill = "white"),
            plot.background = element_rect(fill = "white"),
            panel.border = element_rect(fill = "transparent"),
            axis.ticks = element_blank(),
            plot.margin = unit(c(0.2, 0.02, 0.2, 0.2),"cm")),
    
    # Tax evasion: Wave 3
    ggplot(filter(barplot.waves, Wave == "Wave 3" & DV == "Tax evasion"), 
           aes(x = Answer, y = Share)) + 
      # Grid
      geom_hline(yintercept = seq(5, 45, 5), color = "grey90") +
      
      # Bars
      geom_bar(stat = "identity",
               position = position_dodge(),
               color = "black",
               alpha = 0.66) +
      geom_text(aes(label = paste0(Share, "%\nN=", n)), 
                color = "black", 
                size = 9/.pt, 
                position = position_dodge(0.9), 
                alpha = 0.85, 
                vjust = -0.1, 
                lineheight = 0.8) +
      
      # Right panel
      geom_rect(xmin = 4.65, xmax = 5.15, ymin = 0, ymax = 50, colour = "black", fill = "grey90") +
      annotate("text", x = 4.9, y = 25, label = "Tax evasion", angle = 270, size = 12/.pt) +
      
      # Specs
      scale_y_continuous(breaks = seq(0, 50, 10),
                         labels = c("0%", "10%", "20%", "30%", "40%", "50%")) +
      scale_x_continuous(breaks = c(1, 2, 3, 4), 
                         labels = c("Very\nuncommon", "Rather\nuncommon", "Rather\ncommon", "Very\ncommon")) +
      labs(x = NULL, y = NULL, fill = NULL) +
      coord_cartesian(xlim = c(0.35, 5.15), ylim = c(0, 50), expand = FALSE, clip = "off") +
      geom_hline(yintercept = c(0, 50)) + geom_vline(xintercept = 0.35) +
      theme(text = element_text(color = "black"),
            axis.text.y = element_blank(),
            legend.spacing.y = unit(0.25, "cm"),
            legend.key.size = unit(0.65, "cm"),
            panel.grid = element_blank(), 
            panel.background = element_rect(fill = "white"),
            plot.background = element_rect(fill = "white"),
            panel.border = element_rect(color = "transparent", fill = "transparent"),
            axis.ticks = element_blank(),
            plot.margin = unit(c(0.2, 0.02, 0.2, 0.2),"cm")),
    
    # Shared figure specs
    ncol = 3, widths = c(4.95, 4.95, 5.15),
    nrow = 2)


# Saving the output
ggsave(barplot.waves.fig, filename = "04-figures and models/Fig A1 - barplots per waves.png",
       width = 25, height = 15, units = "cm", dpi = 600)



# Interim cleanup
rm(barplot.waves, barplot.waves.fig)




# _ ROBUSTNESS TESTs ------------------------------------------------------


# __ ALTERNATIVE MEASURE of system abuse ----------------------------------
# NOTE: This robustness test replicates the dependent variable. It uses
#       question 'Many unemployed people do not really try to find a job'
#       (q27_8). This variable is however included only in the first
#       wave of SUPPA panel survey--therefore, only cross-sectional
#       analysis is possible.



# ___ Table S5 ------------------------------------------------------------

ols.unp.benefit <- feols(q27_8 ~ benefit +
                           sw0(male + alder + Q34 + as.factor(persinnt12) + as.factor(utdanning)),
                         weights = ~vekt, 
                         data = na.omit(SUPPA[c("q27_8", 
                                                "benefit", 
                                                "male", 
                                                "alder", 
                                                "Q34", 
                                                "persinnt12", 
                                                "utdanning", 
                                                "id", 
                                                "vekt")]))


# Summary
etable(ols.unp.benefit)



# Exporting the table
write_tableHTML(
  tableHTML(etable(ols.unp.benefit, tex = FALSE, digits = 3)), 
  file = "04-figures and models/Tab S5 - alt unemp indicator - between-group differences.html")



# Interim cleanup
rm(ols.unp.benefit)





# _ FLIPPED MODELS --------------------------------------------------------

# __ Table S6: Cross-section (flipped DVs) --------------------------------

# Models S3 & S4: Tax evasion
ols.flip.tax <- feols(Q30b_2 ~ benefit +
                        sw0(
                          as.factor(Periode) + 
                            male + alder + Q34 + as.factor(persinnt12) + as.factor(utdanning)
                        ),
                      cluster = ~id, 
                      weights = ~vekt, 
                      data = na.omit(SUPPA[c("Q30b_2",
                                             "benefit", 
                                             "Periode", 
                                             "male", 
                                             "alder", 
                                             "Q34", 
                                             "persinnt12", 
                                             "utdanning", 
                                             "id", 
                                             "vekt")]))


# Models S5 & S6: Benefit abuse
ols.flip.benefit <- feols(Q30b_1 ~ inc_increase + inc_decrease +
                            sw0(
                              as.factor(Periode) + 
                                male + alder + Q34 + as.factor(persinnt12) + as.factor(utdanning)
                            ),
                          cluster = ~id, 
                          weights = ~vekt, 
                          data = na.omit(SUPPA[c("Q30b_1", 
                                                 "inc_increase", 
                                                 "inc_decrease", 
                                                 "Periode", 
                                                 "male", 
                                                 "alder", 
                                                 "Q34", 
                                                 "persinnt12", 
                                                 "utdanning", 
                                                 "id", 
                                                 "vekt")]))


# Summary
etable(ols.flip.tax, ols.flip.benefit)



# Exporting the table
write_tableHTML(
  tableHTML(etable(ols.flip.tax, ols.flip.benefit, tex = FALSE, digits = 3)), 
  file = "04-figures and models/Tab S6 - flipped DVs - between-group differences.html")





# ___ Table S7: Within-individual (flipped DVs) ---------------------------

# TAX EVASION: Models
# BASE model
fix.tax.flip.base <- feols(Q30b_2 ~ benefit | id + Periode,  
                           cluster = ~id, 
                           weights = ~vekt, 
                           data = SUPPA)

# Effect heterogeneity
fix.tax.flip.inc_grp <- feols(Q30b_2 ~ benefit * income_grp | id + Periode,  
                              cluster = ~id, 
                              weights = ~vekt, 
                              data = SUPPA)




# BENEFIT ABUSE: Models
# BASE model
fix.benefit.flip.base <- feols(Q30b_1 ~ persinnt12 | id + Periode,  
                               cluster = ~id, 
                               weights = ~vekt, 
                               data = SUPPA)

# Effect heterogeneity
fix.benefit.flip.inc_grp <- feols(Q30b_1 ~ persinnt12 * income_grp | id + Periode,  
                                  cluster = ~id, 
                                  weights = ~vekt, 
                                  data = SUPPA)



# Output
etable(fix.tax.flip.base, fix.tax.flip.inc_grp, fix.benefit.flip.base, fix.benefit.flip.inc_grp)




# Exporting the table
write_tableHTML(
  tableHTML(etable(fix.tax.flip.base, fix.tax.flip.inc_grp, fix.benefit.flip.base, fix.benefit.flip.inc_grp, tex = FALSE, digits = 3)), 
  file = "04-figures and models/Tab S7 - flipped DVs - within-group.html")




# Final cleanup
rm(ols.flip.tax, ols.flip.benefit)
rm(fix.tax.flip.base, fix.tax.flip.inc_grp, fix.benefit.flip.base, fix.benefit.flip.inc_grp)
rm(SUPPA)


### END